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INTRODUCTION 

Computation of the transonic flow around airfoils and wings 
using finite-difference and finite-volume methods requires 
fine grids and large computational domains around the source 
of disturbance. The outer boundary of the computational 
domain is usually placed at several chord lengths away from 
the inner boundary. Moreover, special treatment is required 
at the outer boundary to approximately satisfy the farfield 
boundary conditions. Using different levels of 

approximations, lnvlscid computational schemes have been 

developed based on the Transonic Small Perturbation (TSP) 
equation (Murman and Cole* ; Edwards, Bland and Siedel^), Full 
Potential (FP) equation (Steger and Lomax' 3 ; Garabedlan and 

Korn 4 ; Jameson^) and Euler equation (Jameson**). These schemes 
require large capacity of computer memory to handle the large 
number of grid points and the associated flow variables. They 
also require large CPU time to obtain converged solutions due 

to the large number of iteration cycles usually it is of 

order a thousand. 

The potential equations can be used for flows with weak 
schocks since the entropy increase and vorticity production 
across these shocks are small. For strong shocks, the 
irrotational and isen tropic flow assumptions are invalid and 
one cannot use the potential flow equations across or behind 
the shocks. For these flows, one has to correct the potential 
equations in order to include the entropy jump across the 
shock wave. Methods of this type have been developed by Hafez 
and Lovell^, Fuglsang and Williams® and Whitlow, Hafez and 
Osher*. Alternatively, one has to use the computationally 


more expensive Euler equations. 

The Integral Equation Formulation (IEF) using the TSP 
(Piers and Sloof 10 ; Tseng and Morino 11 ) or the FP (Kandil and 
Yates^; Oskam*^; Erickson and Strande 1 ^; Sinclair*^; Kandil 
and Hu***) equations represents an alternative to the finite 
difference and finite volume methods to treat transonic flows 
with weak shocks. With the IEF, the farfield boundary condi- 
tions are automatically satisfied and only a small computa- 
tional domain is needed around the source of disturbance. 
Moreover, the accuracy of the method depends on the evaluation 
of integrals rather than derivatives and hence coarse grids 
can be used within the small computational domain. Because of 
these obvious advantages of the methods which are based on the 
IEF, it is highly desirable to fully develop these methods and 
extend them to treat transonic flows over a wide range of Mach 
numbers. 

In this paper, we present a transonic integral equation 
method which is based on the full potential equation, and 
couple the method with embedded Euler domains to treat strong 
shocks. The method is extensively applied to different 
airfoil sections and Mach numbers. The results are compared 
with experimental data and other results which were obtained 
by using finite-difference and finite-volume methods with TSP, 
FP and Euler equations. 

FORMULATION 


Full Potential Equation 

For transonic flows with weak shocks, the dimensionless 
governing equations of the two-dimensional, steady, potential 
flow are given by 
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G - T (p * ®* + p y V <2> 

P - [1 +^M 2 (1 -®y)] l/rl . (3) 


where <5 is the total velocity potential, p the density, 
the frees tream Mach number, y the ratio of specific heats and 
the subscripts x and y refer to the derivatives. It should be 
noted that G represents the total compressibility in the flow; 
The characteristic parameters are the frees tream velocity 
(U ) and density (p ) and the airfoil chord length (A). 
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The boundary conditions on the airfoil g, away from g and at 
the trailing edge TE are given by 
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Vg/|vg|, e,, Is 

a unit vector parallel to U 

and 


ACp is the pressure jump. The pressure coefficent is given by 
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The formal solution of Eq. (1) in terms of the velocity field 
V® with explicit contribution of the shock surface S is 
obtained as 


v®<*,y) - 1. + ^ i V »> ~ I 7 ' 11 !; 3 ds 

g 6 (x-C) + (y— n) 



where q g and y are the airfoil surface distributions of 
sources and vorticity; respectively, and q^ is the source 
strength of the shock surface. In the shock-capturing scheme, 
the last integral term In Eq. (8) is dropped, since this term 
is Included in the third integral term of the equation. In 
the shock capturing- shock fitting scheme, the last integral 
term corresponding to the shock .surface is retained when shock 
fitting is used. For shock fitting, the following equations 
are used to determine the shock strength q^, properties behind 
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the shock p^, V 2tx* V 2t an< * t * ie or * entat * on the segments 
forming the shock surface 0: 


2 V 


«1. - V 2«> 


In 


Y+l 


(1 - 


M 


In 


2n 


(Y-D Mj n + 2 

(T+l) «L 


In 


M ln > 1 


(9) 


( 10 ) 


V a V 
v 2t It 


( 11 ) 


( T*> M in 

(Y-D M^ n + 2 


( 12 ) 


p » sin [ 


-1 rl.2 sing slnQ . 1 


cos (p-9) 


-41 


1/2 


m: 


(13) 


M, 


r± ni 

MA 2 • "u ■ 7 ®1 • V P 1 2 


(14) 


where the subscripts 1 and 2 refer to the conditions ahead and 
behind the shock, respectively, while n and t refer to the 
normal and tangential directions to the shock surface S, 
respectively, and 9 is the relative direction of the flow 
behind the shock to that ahead of the shock. 


Euler Equations 

For strong shocks, an embedded computational domain is 
constructed around the shock which has been preliminary found 
by the Integral solution with shock capturing only. With the 
boundary and initial conditions found from the Integral 
solution, the unsteady conservative form of the Euler 
equations are solved in this limited domain with psuedo time 
marching. The dimensionless conservative form of these 
equations is given by 


M + M + 2L m 

dt 5x 5y 


0 


(15) 


where the flow vector field q and the flux components E and 
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F are given by 


q 

E 

F 


[p, pu, pv, pe] t 


(16) 

[pu, pu 2 + p, puv , 

puh] t 

(17) 

[pv, puv, pv 2 + p. 

pvh] C 

(18) 


The total energy and enthalpy per unit mass are given by 

6 “ (y+i) p + + y2 ^ 2 * h " e + p/p (19) 


Since we are interested in the steady flow solution only, the 
energy equation [last elements in Eqs. (16)-(18)], which is a 
differential equation, is replaced by the algebraic steady 
form which states that the total enthalpy is constant. Hence, 
the energy equation is replaced by 

p - (p/y) [1/M* + (y-1) (1 - u 2 - v 2 )/2] (20) 

METHOD OF SOLUTION 

Shock-Capturing Shock-Fitting (SCSF) Scheme 

The basic difference between the incompressible Integral 
Equation Solution and the transonic Integral Equation Solution 
are the additional third and fourth integral terras of Eq. (8). 
In the shock capturing part of the scheme, the fourth integral 
terra is dropped while in the shock fitting part of the scheme 
this term is taken into account. 

The SCSF-scheme is an iterative scheme due to the 
nonlinearity of the third and fourth integral terms. The 
iterative scheme is described below: 

Neglecting the fourth integral term and setting G ■ 0, a 
standard panel computation is used to obtain q g and/or y . 
The q and y are piecewise linear distributions on each 
surface panel ^and they are defined in terms of their nodal 
values. Initial values of G are calculated at the centroids 
of the field elements by using the linear compressibility 

relation G - M 2 u . where u„ Is the x-derivative of the x- 

CD X X 

component of the velocity. The centroidal value G represents 
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the G value for the field element. Equation (8) Is then used 
to enforce the surface boundary conditions, Eqs. (5) and (7) 
to find new q g and/or y g . The density p and nonlinear 
compressibility G are calculated by using Eqs. (3) and (2). A 
type finite-difference expression is used to calculate p and 
p depending on the type of the centroidal point-subsonic or 
supersonic. Once the G values are obtained, the surface 
boundary conditions are satisfied again. The iterative 
procedure is continued until the shock location is fixed. 
This is the shock capturing part of the scheme. 

Shock panels are then introduced at the shock location, 
the fourth itegral term of Eq. (8) is now taken into account, 
Eqs. (9) and (13) are used to calculate q 3 and p and Eqs. 

(10)-(12) are used to cross the shcok panel. The iterative 
procedure is continued as before with the exception of dealing 
with the shock panels as explained. Convergence is achieved 
once the surface pressure converges. This is the shock 
fitting part of the scheme. 

Integral Equation With Embedded Euler (IEEE) Scheme 
In this scheme, the shock capturing part of the SCSF-scheme is 
used to locate the shock. Once the shock is captured, a fine 
grid is constructed within a small computational region around 
the shock where a finite-volume Euler scheme is used. The 
basic f ini te-volume equation is obtained by integrating Eq. 
(15) over x and y to obtain 


// ^ dA + <J) (E dy + F dx) - 0 (21) 

Equation (21) is then applied to each cell of the embedded 
grid of the Euler domain. The resulting difference equation 
is 

4 

(§J) A A + Z (E A y + F A x ) ■ 0 (22) 

5t u 13 t-i r r 

where AA is the cell area, r refers to the cell-side number 
and the integer subscripts i, j refer to the centroidal 
values. The Euler solver is a central-difference finite- 
volume method which uses four-stage Runge-Kutta time stepping 
with explicit second- and fourth-order dissipation terms. The 
details of this solver are given in reference [17]. 

The boundary and initial conditions for the Euler domain 
are obtained from the Integral equation solution which is 
interpolated on the Euler domain grid. The Euler solver is 
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then used to capture the shock and calculate the flow vector 
field q. It should be emphasized here that the downstream 
boundary condition must be_ updated while Euler calculations 
are excuted. Fixing the q values of the Euler domain, the 
IE calculation is used to update the boundary conditions for 
the Euler domain. The iterative procedure is repeated until 
convergence is achieved. 

COMPUTATIONAL EXAMPLES 

In this section, we present applications of the SCSF- and 
IEEE-scheraes to the NACA 0012 and 64A010A. According to the 
convergence study using different sizes of the IE 
computational domain, which was presented by the authors 

(Kandll and Hu 1 ^), a computational domain of 2 x 1.5, in the x 
and y directions, has been used around the airfoil in all the 
following applications. A rectangular grid of 64 x 60 has 
been used for the IE computation. The grid is clustered in 
the leading edge, the shock region and near the airfoil 
surface. Since the third Integral terra of Eq. (8) is 
computationally expensive, its computation with constant G 
distribution has been restricted to the nearfield computation. 
For the farfield computation, this term is replaced by the 
equivalent lumped source term at its centroid. With 
sufficient accuracy, it has been computationally determined 
that the nearfield distance from. the centroid is < 0.5. 

Figure 1 shows the results of the SCSF-scheme for NACA 
0012, M * 0.8 and a ■ 0° , along with comparisons with the 

IQ 

computational results of Garabedian, Korn and Jameson , and 
the experimental data taken from reference 19. The SCSF 
scheme took 12 iteration cycles of shock capturing (SC) and 13 
cycles of shock fitting (SF) to achieve convergence. 

Figure 2 shows the results of the IEEE-scheme for the 
same case along with a comparison with the computational 
results of Jameson** who also used the finite-volume Euler 
scheme with four-stage Runge Kutta time stepping. In the 
present IEEE-scheme, the embedded Euler domain has a size of 
0.5 x Q.6 around the shock region with a grid of 25 x 30. 
This case took 10 iteration cycles of SC, 250 time cycles of 
Euler iterations to achieve a residual error of 10~^ and 5 IE 
cycles to update the Euler domain boundary conditions. 

Figures 3 and 4 show the results of the SCSF- and IEEE- 
schemes for NACA 64A010A, M o a 0.796, a ** 0° along with 
comparisons with the computational results of Edwards, Bland 
and Seidel 2 who used the TSP equation, and the experimental 
data taken from reference 2.® With the SCSF-scheme, the 
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numbers of SC and SF iteration cycles to achieve convergence 
are the same as those of the case presented in Fig. 1. With 
the IEEE-scheme, the embedded Euler domain has a size of 0.7 x 
0.6 with a grid size of 35 x 30. This case. Fig. 4, took 10 
iteration cycles of SC, 130 time cycles of Euler iterations to 
achieve a residual error of 10"^ and 3 IE cycles to update the 
Euler domain boundary conditions. 

Figures 5 and 6 show the results of the SCSF- and IEEE- 
schemes for the lifting case of NACA 0012, « 0.75 and 

a » 2° along with the computational results of Steger and 
Loraax^ , and the experimental data taken from the same 
reference. The size of the grids and the number of iteration 
cycles used to achieve convergence are the same as those of 
the cases given in Fig. 1 and 2. 

Figure 7 shows the results of the IEEE for NACA 0012, 
■ 0.812 and cc ■ 0° along with the experimental data of 
reference 18. In Fig. 8, the results of the IEEE for NACA 
0012, ■ 0.82 and a ■ 0° are shown along with the three- 

dimensional solution at the wing root chord of Tseng and 

Morino^, who use the IE for the TSP, and the experimental 
results of reference 20. The size of the embedded Euler 

domain for these cases is 0.8 x 0.8 end the computational grid 
is 40 x 40. 

CONCLUDING REMARKS 

Two transonic computational schemes which are based on the 
Integral Equation Formulation of the full potential equation 
have been presented. The first scheme is a Shock Capturing- 
Shock Fitting (SCSF) scheme which uses the full potential 

equation throughout with the exception of the shock wave where 
the Rankine-Hugoniot relations are used to cross and fit the 
shock. The second scheme is .an Integral Equation with 

Embedded Euler (IEEE) scheme which uses the full potential 

equation with an embedded region where Euler equations are 
used. The two schemes are applied to several transonic 

airfoil flows and the results have been compared with numerous 
computational results and experimental data. The two schemes 
are nevertheless efficient as compared to the other existing 
schemes which use finite-difference or finite-volume methods 
throughout large computational domains with fine grids. The 
SCSF-scheme is restricted to flows with weak shocks, while the 
IEEE-scheme can handle strong shocks. Currently, the IEEE 
scheme is applied to other transonic flows with strong shocks 
as well as to unsteady pitching oscillations. 
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